Counting the socks¶

Introduction¶

Came across an amusing little probability problem the other day, from Fifty Challenging Problems in Probability with Solutions, by Frederick Mostell, with connections to some deep mathematics

The question is: given you have red and black socks in a drawer, how many of each do you need, to have a 50% chance of getting two red socks when picking at random from the drawer (sampling without replacement)r

In [1]:

import math
import numpy as np
import sympy

Analysis¶

Given we have $\require{cancel}$ $r$ red socks and $b$ black socks, with $n$ socks in total, the probability of getting two red socks is:

$$\frac{r * (r-1)} { n * (n-1)}$$

so we want integers $r$ and $n$ such that

$$\frac{r * (r-1)} { n * (n-1)} = \frac{1}{2}$$

Re-arranging terms, we get

$$ 2*r^2 -2*r = n*(n-1)$$

If $A*x^2+B*x+C = 0$, we have solutions $\frac{-B \pm \sqrt{B^2-4*A*C}}{2*A}$

Here, we have $A=2, B=2, C=-n*(n-1)$, so

$$ r = \frac{2 \pm \sqrt{4+4*2*n(n-1)}}{2*2}$$ or $$ r = \frac{2 \pm \sqrt{4+8*n(n-1)}}{2*2}$$

so for $r$ to be an integer, we must have $[1+2*n*(n-1)] $ to be a square

Mathematical connections¶

This is an example of a Diophantine equation:

a Diophantine equation is a polynomial equation with integer coefficients, for which only integer solutions are of interest.

Diophantine equations have been studied for very long time: Diophantus

may have lived anywhere between 170 BCE, roughly contemporaneous with Hypsicles, the latest author he quotes from, and 350 CE, when Theon of Alexandria quotes from him

Pierre de Fermat famously stated his "Last Theorem" around 1637 in the margin of a copy of the ancient Greek text Arithmetica by Diophantus. This problem is also connected to Pell's equation: from OEIS-

... this is equivalent to the Pell equation A(n)^2-2*B(n)^2 = -1

Exploration¶

First I search out to find the first few candidate values, where $[1+2*n*(n-1)] $ is a a square

In [2]:

for i in range(2,5000):
    x = math.sqrt(1+2*i*(i-1))
    
    # test for integer solution
    if( int(x) == x):
        red = int((1+int(x))/2)
        print(f'Sock count {i:4}, with {red:3} red socks')
    #end if
#end for

Sock count    4, with   3 red socks
Sock count   21, with  15 red socks
Sock count  120, with  85 red socks
Sock count  697, with 493 red socks
Sock count 4060, with 2871 red socks

Checking, we get the probability of getting 2 red socks from 3 red socks out of 4 in the drawer $$ \frac{\cancel{3}}{4}*\frac{2}{\cancel{3}} = \frac{1}{2}$$

Checking, we get the probability of getting 2 red socks from $15 ( = 3*5)$ red socks out of $21 ( = 3*7)$ in the drawer

$$ \frac{\cancel{3}*\cancel{5}}{\cancel{3}*\cancel{7}}*\frac{2*\cancel{7}}{4*\cancel{5}} = \frac{1}{2}$$

Now it turns out that $4, 21, 120, 697, ...$ is a well known sequence, listed in The On-Line Encyclopedia of Integer Sequences [https://oeis.org/A046090], as a solution to this very problem.

Recurrance relations¶

The OEIS gives a recurrance relation for computing the next valid sock count (given the two previous sock count), but I thought I would find it for myself. Let us assume that any valid value depends only on the previous two values (a la Fibonacci) and a constant. If $a_{n}$ is the n-th value in the series:

$$ a_{n} = A*a_{n-1} + B*a_{n-2} + C$$

so we have

$$ \begin{align*}4060 &= A*697 + B*120 + C\\ 697 &= A*120 + B*21 + C\\ 120 &= A*21 + B*4 + C\\ \end{align*} $$

To solve for $A, B, C$, we use numpy linear algebra routines

In [3]:

b  = [4060, 697, 120]

matrix = np.array([ [697,120,1], [120,21,1], [21,4,1] ])

sol = np.linalg.solve(matrix,b)

print(f'Recurrance relation: a(n) = {round(sol[0]):+}*a(n-1)   {round(sol[1]):+}*a(n-2)  {round(sol[2]):+}')

Recurrance relation: a(n) = +6*a(n-1)   -1*a(n-2)  -2

So to get the first few values of valid sock count values

In [4]:

# initial value in series
a_n_1 = 4
a_n = 21

values = [4,21]

# find next 30 values
for i in range(30):
    # compute next value in series
    x = 6*a_n - a_n_1 -2
    # add it to out list
    values.append(x)
    # update recurrance values
    a_n, a_n_1 =   x,a_n
#end for

# display values nicely
for i,v in enumerate(values):
    print(f' {i+1:2}:  {v:40,}')
#end for

  1:                                         4
  2:                                        21
  3:                                       120
  4:                                       697
  5:                                     4,060
  6:                                    23,661
  7:                                   137,904
  8:                                   803,761
  9:                                 4,684,660
 10:                                27,304,197
 11:                               159,140,520
 12:                               927,538,921
 13:                             5,406,093,004
 14:                            31,509,019,101
 15:                           183,648,021,600
 16:                         1,070,379,110,497
 17:                         6,238,626,641,380
 18:                        36,361,380,737,781
 19:                       211,929,657,785,304
 20:                     1,235,216,565,974,041
 21:                     7,199,369,738,058,940
 22:                    41,961,001,862,379,597
 23:                   244,566,641,436,218,640
 24:                 1,425,438,846,754,932,241
 25:                 8,308,066,439,093,374,804
 26:                48,422,959,787,805,316,581
 27:               282,229,692,287,738,524,680
 28:             1,644,955,193,938,625,831,497
 29:             9,587,501,471,344,016,464,300
 30:            55,880,053,634,125,472,954,301
 31:           325,692,820,333,408,821,261,504
 32:         1,898,276,868,366,327,454,614,721

We can see that we get some very large integers very quickly! In manipulating these numbers (especially taking square roots), we move to a representation as powers of prime factors.

Checking a large value¶

We are dealing with very large integers, so it is more convenient to deal with numbers encoded as powers of their prime factors. We define some helper functions.

We use sympy to factor our numbers

In [5]:

def sqrt_factors(d):
    '''
    sqrt_factors: given a dictionary of prime factors of a number, return the prime factors of the sqrt
    '''
    d2 = {}
    for k in d.keys():
        d2[k] = round(d[k]/2)
    #end for
    return d2
#end sqrt_factors

def sqrt_value(d):
    '''
    sqrt_value: given a dictionary of prime factor of a number, return the sqrt value
    '''
    d2 = sqrt_factors(d)
    value = 1
    for k in d2.keys():
        for p in range(d2[k]):
            value = value * k
        #end for
    #end for
    return value
#end sqrt_value

def div_by_4(n):
    '''
    div_by_4: given a number, return (number / 4 )

    assumes number is divisible by 4, so reduce power of 2 in factor list by 2
    '''    
    d = sympy.factorint(n)
    d[2] = d[2]-2
    
    # multiply out factors
    value = 1
    for k in d.keys():
        value = value * k**d[k]
    #end for
    return value
#end div_by_4

In [6]:

def show_factors(factors):
    '''
    show_factors: give a dictionary of prime factors, return a string depiecting products of powers
    '''
    s = ""
    first = True
    for k in factors.keys():
        # * only appears between factors
        if(first):
            s = s+ f'{k}^{factors[k]}'
            first = False
        else:
            s = s+ f' * {k}^{factors[k]}'
        #end if
    #end for
    return s
#end show_factors

In [7]:

def get_r_value(n):
    '''
    get_r_value: given the total sock count, find red sock count
    '''
    x = 4 + 8*(n*n-n)
    y = sqrt_value(sympy.factorint(x))
    r = div_by_4(2+y)
    return r
#end get_r_value

In [8]:

def show_prob_calculation(v):   
    '''
    show_prob_calculation: given a valid sock count, show the factors of (r/n)*((r-1)/(n-1))
    '''
    
    z = v
    # get red sock count
    r = get_r_value(z)
    print ( f'Red sock count {r:,}, Total sock count {z:,}')

    # show  (r/n)*((r-1)/(n-1)) as product of powers of prime factores
    print(
        '(' + show_factors(sympy.factorint(r-1)) + ') * (' +
        show_factors(sympy.factorint(r)) + ')'
         
         )
    print('-'*80)
    print(
        '(' + show_factors(sympy.factorint(z-1)) + ') * (' +
        show_factors(sympy.factorint(z)) + ')'
         
         )
    print('\n\n')
    #end for
#end show_prob_calculation

Now we can show how the factors of $(r/n)/((r-1)/(n-1))$ cancel out, leaving $\frac{1}{2}$, for our largest computed valid total sock count

In [9]:

show_prob_calculation(values[-1])

Red sock count 1,342,284,446,191,393,385,645,665, Total sock count 1,898,276,868,366,327,454,614,721
(2^5 * 3^1 * 7^1 * 17^1 * 23^1 * 353^1 * 577^1 * 665857^1 * 37667521^1) * (5^1 * 257^1 * 1409^1 * 5741^1 * 2448769^1 * 52734529^1)
--------------------------------------------------------------------------------
(2^6 * 3^1 * 5^1 * 17^1 * 577^1 * 5741^1 * 665857^1 * 52734529^1) * (7^1 * 23^1 * 257^1 * 353^1 * 1409^1 * 2448769^1 * 37667521^1)

Conclusion¶

I wouldn't have guessed that such a simpe-sounding problem would have such deep mathematics roots

There is a saying (when you want to decline an invitation) "That's the day I'm re-organizing my sock drawer": based upon the above, this could indeed take some time!

Reproducability information¶

In [10]:

%load_ext watermark

In [11]:

%watermark

Last updated: 2026-08-20T19:55:35.266650+10:00

Python implementation: CPython
Python version       : 3.11.7
IPython version      : 8.20.0

Compiler    : MSC v.1916 64 bit (AMD64)
OS          : Windows
Release     : 10
Machine     : AMD64
Processor   : Intel64 Family 6 Model 170 Stepping 4, GenuineIntel
CPU cores   : 22
Architecture: 64bit

In [12]:

%watermark -h -iv -co

conda environment: base

Hostname: INSPIRON16

numpy     : 1.26.4
json      : 2.0.9
pandas    : 2.1.4
ipywidgets: 7.6.5
xarray    : 2025.4.0
sympy     : 1.12
sys       : 3.11.7 | packaged by Anaconda, Inc. | (main, Dec 15 2023, 18:05:47) [MSC v.1916 64 bit (AMD64)]

In [ ]:

Comment

Counting the socks